/* CenterMass.h	-- Massimo Marchi
 *
 * $Header$
 * $Log$
 */
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>
#include <vector>
#include <algorithm>

#include "MyUtilClass.h"
#include "Quaternions.hpp"
#include "WignerDMatrices.hpp"
#include <sstream>
#ifndef _CENTERMASS_H
#define _CENTERMASS_H

using namespace std;
using namespace DVECT;
using namespace MATRIX;
using vectorvd=vector<vector<Dvect> >;
using vectorvi=vector<vector<int> >;
using Quaternions::Quaternion;
using SphericalFunctions::WignerDMatrix;
using  vvvector_d=vector<vector<vector<Dvect> > >;
using  vvector_d=vector<vector<Dvect> >;
using  vvector_i=vector<vector<int> >;



class CenterMass{
protected:
  vector<Dvect> cm;
  vector<vector<Dvect> > cmm;
  Dvect Axis;
  Matrix CO,OC;
  static double time_c;
  vector<Quaternion> Q_t;
  virtual void WriteIt(ostream &);
  virtual void ReadIt(ifstream &);
  bool DoneOnce;
  Matrix MatrixfromQ(Quaternions::Quaternion );
  void RotateBack(vector<vector<Dvect> > & );
public:
  CenterMass(): DoneOnce{false} {cm.clear();cmm.clear();};
  static void setTime(double tt){time_c=tt;};
  void setCOs(Matrix co, Matrix oc){ CO(co);OC(oc);}
  void setAxis(double x, double y, double z){Axis[XX]=x;Axis[YY]=y;Axis[ZZ]=z;}
  void setup(size_t m){
	  cm=vector<Dvect>(m,{0.0,0.0,0.0});
	  cmm=vectorvd(m);
	  Dvect t{0.0,0.0,0.0};
	  vector<Dvect> v=vector<Dvect>(8,Dvect{0.0,0.0,0.0});
	  Q_t=vector<Quaternion>(m,Quaternion{1.0,0.0,0.0,0.0});
  }
  Dvect & getAxis(){return Axis;}
  static double getTime(){return time_c;};
  virtual void operator()(vector<Dvect> &, vector<vector<Dvect> >&,  vector<Quaternion> &);
  virtual void operator()(vector<Dvect> &, vector<Quaternion> &);
  virtual void operator()(vector<Dvect> &);
  virtual void Diffusion(ofstream & fout){};

  size_t Size(){return cm.size();};
  Dvect & operator[](const int i) {return cm[i];};
  vector<Quaternion> & getQ(){return Q_t;}
  vector<Dvect> & getR(){return cm;};
  virtual ~CenterMass(){};
  virtual void RefCoord(rvec * y, vector<string> &y0,vector<vector<int> > & z, vvector_i & z0){};

  friend ostream & operator<<(ostream &, CenterMass &);
  friend ifstream & operator>>(ifstream &, CenterMass & );
};


#endif
